# Appendix figures for: 
# Class, Policy Attitudes and U.S. Presidential Voting in the Post-Industrial Era: 
#  The Importance of Issue Salience
# Franko & Witko
# Political Research Quarterly
# Date: 05/25/22


# Install and load required packages.
if (!require(tidyverse)) install.packages("tidyverse")
library(tidyverse)
if (!require(haven)) install.packages("haven")
library(haven)
if (!require(janitor)) install.packages("janitor")
library(janitor)
if (!require(gridExtra)) install.packages("gridExtra")
library(gridExtra)


### Change path to directory where replication files are stored. ###
setwd("~/Dropbox/Class/ClassPolicyAttitudesVotingPaper/prq-final-data")


##################################################
# Set theme options.
##################################################

color.background <- "#F0F0F0"
color.grid.major <- "#DBDBDB"
fte_theme <- function() {
  theme(plot.background = element_rect(fill = color.background)) + 
    theme(panel.background = element_rect(fill = color.background, 
                                          colour = color.background)) +
    theme(panel.grid.major=element_line(color = color.grid.major, size = .25)) +
    theme(panel.grid.minor = element_blank()) +
    theme(axis.ticks = element_blank()) +
    theme(axis.title.x = element_text(vjust = -0.5), 
          axis.title.y = element_text(vjust = 1.2)) +
    theme(plot.title = element_text(hjust = 0.5))
}


##################################################
# Descriptive plot.
##################################################

# Load data.
gss <- read_dta("GSS_Class_wEGP_clean.dta")

gss.0616 <- filter(gss, year>=2006, !is.na(egp5))
count(gss.0616, egp5)

# Labels.
egp5_labs <- c("1"="Upper Class", "2"="Upper Middle Class", "3"="Middle-Class Serv.", 
               "4"="Middle-Class Man.", "5"="Working Class")

# Income and education within EGP class categories.
figEGP.inc <- ggplot(filter(gss.0616, !is.na(egp5), !is.na(inc5)), 
                     aes(factor(inc5), weight = wtss)) +
  geom_histogram(aes(y = ..count.. / tapply(..count..,..PANEL.., sum)[..PANEL..]), 
                 stat = "count", fill = "#66CC99") +
  labs(y = "", x = "Respondent Income") +
  facet_wrap( ~ egp5, ncol = 2, labeller = labeller(egp5 = egp5_labs)) +
  fte_theme() 


figEGP.edu <- ggplot(filter(gss.0616, !is.na(egp5), !is.na(edu5)), 
                     aes(factor(edu5), weight = wtss)) +
  geom_histogram(aes(y = ..count.. / tapply(..count..,..PANEL.., sum)[..PANEL..]), 
                 stat = "count", fill = "#619cff") +
  labs(y = "", x = "Respondent Education") +
  facet_wrap( ~ egp5, ncol = 2, labeller = labeller(egp5 = egp5_labs)) +
  fte_theme() 


# Combine plots.
pdf("figB1.pdf", width=6, height=3)
grid.arrange(figEGP.inc, figEGP.edu, ncol=2)
dev.off()



##################################################
# Alternative models excluding party ID and 
# ideology as covariates.
##################################################


# Labels.

class_labs <- c("edu5"="Edu. low -\nedu. high", 
                "egp1.5"="Working\nclass -\nupper class", 
                "egp2.5"="Working\nclass -\nupper middle", 
                "inc5"="Inc. low -\ninc. high", 
                "subclass"=
                  "Subjective\nlower class -\nupper class", 
                "white"="Nonwhite -\nwhite")

class_labs2 <- c("edu5"="Edu. low -\nedu. high", 
                 "egp1.5"="Working class -\nupper class", 
                 "egp2.5"="Working class -\nupper middle", 
                 "inc5"="Inc. low -\ninc. high", 
                 "subclass"=
                   "Subj. low class -\nupper class", 
                 "white"="Nonwhite -\nwhite")


# Class and voting.

# Low/high difference effects.

# Load data.
vdiffs1 <- read_dta("results/dv_diffs_gss_ALT.dta")
vdiffs2 <- read_dta("results/dv_diffs_anes_ALT.dta")

# Combine difference estimates for plotting.
# Create DV indicator.
vdiffs1$dv <- "Vote Dem. (GSS)"
vdiffs2$dv <- "Vote Dem. (ANES)"

# Combine results.
vdiffs.all <- bind_rows(vdiffs1, vdiffs2)
# Reorder factors.
vdiffs.all$sample <- factor(vdiffs.all$sample, 
                            levels = c("White", 
                                       "Nonwhite", 
                                       "All"))
vdiffs.all$var <- factor(vdiffs.all$var, 
                         levels = c("white", "edu5", 
                                    "inc5", "egp2.5", 
                                    "egp1.5", "subclass"))
vdiffs.all$dv <- factor(vdiffs.all$dv, 
                        levels = c("Vote Dem. (GSS)",
                                   "Vote Dem. (ANES)"))


# Plots.

figDiffsVote <- ggplot(vdiffs.all, aes(y = factor(var), 
                                       x = diff, 
                                       group = sample)) +
  geom_vline(xintercept = 0, color = "gray", 
             linetype = "dashed") +
  ggplot2::geom_errorbarh(aes(xmin = ll, xmax = ul), 
                          size = 1.0, color = "gray", 
                          height = 0, 
                          position=ggstance::position_dodgev(height=0.5)) +
  geom_point(aes(color = sample, shape = sample), size = 2, 
             position=ggstance::position_dodgev(height=0.5)) +
  facet_wrap( ~ dv, ncol = 2) +
  labs(x = "Change in Support", 
       y = "Social Class Characteristics") +
  fte_theme() +
  theme(legend.text = element_text(size = 10),
        legend.title = element_blank(),
        legend.key = element_rect(color = color.background), 
        legend.position = "bottom",
        strip.text.x = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm')) +
  scale_y_discrete(labels = class_labs) +
  scale_shape_discrete(breaks = c("All", 
                                  "Nonwhite", "White")) +
  scale_color_discrete(breaks = c("All", 
                                  "Nonwhite", "White")) +
  theme(plot.background = element_rect(fill = "white"))


pdf(file = "figF1.pdf", width=7.0, height=4.5)
figDiffsVote
dev.off()



# Over time effects.
# Voting.

# Load data.
# GSS.
vgss.egp <- read_dta("results/dv_egp1_gss_T_ALT.dta")
vgss.subclass <- read_dta("results/dv_subclass_gss_T_ALT.dta")
vgss.inc <- read_dta("results/dv_inc5_gss_T_ALT.dta")
vgss.edu <- read_dta("results/dv_edu5_gss_T_ALT.dta")

# ANES
vanes.egp <- read_dta("results/dv_egp1_anes_T_ALT.dta")
vanes.subclass <- read_dta("results/dv_subclass_anes_T_ALT.dta")
vanes.inc <- read_dta("results/dv_inc5_anes_T_ALT.dta")
vanes.edu <- read_dta("results/dv_edu5_anes_T_ALT.dta")

# Create DV indicator.
vgss.egp$dv <- "Vote Dem. (GSS)"
vanes.egp$dv <- "Vote Dem. (ANES)"

# Combine estimates for plotting.
# First need to rename and create variables.
vegplst <- list(vgss.egp, vanes.egp)
vegplst <- vegplst %>%
  map(~ rename(., b = b_egp)) %>%
  map(~ rename(., b_hi = b_egp_hi)) %>%
  map(~ rename(., b_lo = b_egp_lo)) %>%
  map(~ mutate(., var = ifelse(egp5 == 1, "Upper Class", 
                               "Upper Mid. Class")))

# Combine list of data frames.
time.vote <- bind_rows(vegplst)

# Create DV indicators and combine other data frames.
vglst <- list(vgss.subclass, vgss.inc, vgss.edu)
vglst <- vglst %>%
  map(~ mutate(., dv = "Vote Dem. (GSS)"))

valst <- list(vanes.subclass, vanes.inc, vanes.edu)
valst <- valst %>%
  map(~ mutate(., dv = "Vote Dem. (ANES)"))

time.vote <- bind_rows(time.vote, vglst, valst)

# Rename variables.
time.vote <- 
  time.vote %>% 
  mutate(var = recode_factor(var, 
                             "Upper Class" = "Upper Class",
                             "Upper Mid. Class" = "Upper Mid. Class",
                             "subclass" = "Subjective Class",
                             "inc5" = "Income",
                             "edu5" = "Education", 
                             .ordered = TRUE))

# Reorder DV.
vdv.ord <- c("Vote Dem. (GSS)",
             "Vote Dem. (ANES)")
time.vote$dv <- factor(time.vote$dv, 
                       levels = vdv.ord)



# Plots.

figTimeVote <- ggplot(time.vote, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "#00ba38", size = 2) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF2.pdf", width=6, height=6)
figTimeVote
dev.off()



# Class and policy.

# Low/high difference effects.

# Load data.
e.diffs1 <- read_dta("results/ineq_diffs_ALT.dta")
e.diffs2 <- read_dta("results/jobs_diffs_ALT.dta")
c.diffs1 <- read_dta("results/cultG_diffs_ALT.dta")
c.diffs2 <- read_dta("results/cultA_diffs_ALT.dta")
r.diffs1 <- read_dta("results/raceG_diffs_ALT.dta")
r.diffs2 <- read_dta("results/raceA_diffs_ALT.dta")

# White only for race estimates.
r.diffs1 <- filter(r.diffs1, sample=="White") 
r.diffs2 <- filter(r.diffs2, sample=="White") 

# Combine difference estimates for plotting.
# Create DV indicator.
e.diffs1$dv <- "Reduce inequality (GSS)"
e.diffs2$dv <- "Guarantee jobs/inc. (ANES)"
c.diffs1$dv <- "Culture policy (GSS)"
c.diffs2$dv <- "Culture policy (ANES)"
r.diffs1$dv <- "Race policy (GSS)"
r.diffs2$dv <- "Race policy (ANES)"

# Combine results.
pdiffs.all <- bind_rows(e.diffs1, e.diffs2, 
                        c.diffs1, c.diffs2,
                        r.diffs1, r.diffs2)
# Reorder factors.
pdiffs.all$sample <- factor(pdiffs.all$sample, 
                            levels = c("White", 
                                       "Nonwhite", 
                                       "All"))
pdiffs.all$var <- factor(pdiffs.all$var, 
                         levels = c("white", "edu5", 
                                    "inc5","egp2.5", 
                                    "egp1.5", "subclass"))
pdiffs.all$dv <- factor(pdiffs.all$dv, 
                        levels = c("Reduce inequality (GSS)",
                                   "Guarantee jobs/inc. (ANES)",
                                   "Culture policy (GSS)",
                                   "Culture policy (ANES)",
                                   "Race policy (GSS)",
                                   "Race policy (ANES)"))



# Plots.
figDiffsPol <- ggplot(pdiffs.all, aes(y = factor(var), 
                                      x = diff, 
                                      group = sample)) +
  geom_vline(xintercept = 0, color = "gray", 
             linetype = "dashed") +
  ggplot2::geom_errorbarh(aes(xmin = ll, xmax = ul), 
                          size = 1.0, color = "gray", 
                          height = 0, 
                          position=ggstance::position_dodgev(height=0.5)) +
  geom_point(aes(color = sample, shape = sample), size = 2, 
             position=ggstance::position_dodgev(height=0.5)) +
  facet_wrap( ~ dv, ncol = 2) +
  labs(x = "Change in Support", 
       y = "Social Class Characteristics") +
  fte_theme() +
  theme(legend.text = element_text(size = 10),
        legend.title = element_blank(),
        legend.key = element_rect(color = color.background), 
        legend.position = "bottom",
        strip.text.x = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm')) +
  scale_y_discrete(labels = class_labs2) +
  scale_shape_discrete(breaks = c("All", "Nonwhite", "White")) +
  scale_color_discrete(breaks = c("All", "Nonwhite", "White")) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF3.pdf", width=6.5, height=7.0)
figDiffsPol
dev.off()



# Over time effects.
# Redistribution.

# Load data.
# Inequality.
ineq.egp <- read_dta("results/ineq_egp5_T_ALT.dta")
ineq.subclass <- read_dta("results/ineq_subclass_T_ALT.dta")
ineq.inc <- read_dta("results/ineq_inc5_T_ALT.dta")
ineq.edu <- read_dta("results/ineq_edu5_T_ALT.dta")

# Jobs and income.
jobs.egp <- read_dta("results/jobs_egp5_T_ALT.dta")
jobs.subclass <- read_dta("results/jobs_subclass_T_ALT.dta")
jobs.inc <- read_dta("results/jobs_inc5_T_ALT.dta")
jobs.edu <- read_dta("results/jobs_edu5_T_ALT.dta")


# Create DV indicator.
ineq.egp$dv <- "Reduce inequality (GSS)"
jobs.egp$dv <- "Guarantee jobs/inc. (ANES)"

# Combine estimates for plotting.
# First need subset of EGP results.
egplst <- list(ineq.egp, jobs.egp)
egplst <- egplst %>%
  map(~ filter(., egp5==1 | egp5==2)) %>%
  map(~ rename(., b = b_egp)) %>%
  map(~ rename(., b_hi = b_egp_hi)) %>%
  map(~ rename(., b_lo = b_egp_lo)) %>%
  map(~ mutate(., var = ifelse(egp5 == 1, "Upper Class", 
                               "Upper Mid. Class")))

# Combine list of data frames.
time.redist <- bind_rows(egplst)

# Create DV indicators and combine other data frames.
ineqlst <- list(ineq.subclass, ineq.inc, ineq.edu)
ineqlst <- ineqlst %>%
  map(~ mutate(., dv = "Reduce inequality (GSS)"))

jobslst <- list(jobs.subclass, jobs.inc, jobs.edu)
jobslst <- jobslst %>%
  map(~ mutate(., dv = "Guarantee jobs/inc. (ANES)"))

time.redist <- bind_rows(time.redist, ineqlst, jobslst)

# Rename variables.
time.redist <- 
  time.redist %>% 
  mutate(var = recode_factor(var, 
                             "Upper Class" = "Upper Class",
                             "Upper Mid. Class" = "Upper Mid. Class",
                             "subclass" = "Subjective Class",
                             "inc5" = "Income",
                             "edu5" = "Education", .ordered = TRUE))

# Reorder DV.
dv.ord <- c("Reduce inequality (GSS)",
            "Guarantee jobs/inc. (ANES)")
time.redist$dv <- factor(time.redist$dv, 
                         levels = dv.ord)


# Plots.

figTimeRedist <- ggplot(time.redist, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "#619cff", size = 2) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF4.pdf", width=6, height=6)
figTimeRedist
dev.off()



# Over time effects.
# Culture.

# Load data.
# GSS.
cultgss.egp <- read_dta("results/cultG_egp5_T_ALT.dta")
cultgss.subclass <- read_dta("results/cultG_subclass_T_ALT.dta")
cultgss.inc <- read_dta("results/cultG_inc5_T_ALT.dta")
cultgss.edu <- read_dta("results/cultG_edu5_T_ALT.dta")

# ANES.
cultanes.egp <- read_dta("results/cultA_egp5_T_ALT.dta")
cultanes.subclass <- read_dta("results/cultA_subclass_T_ALT.dta")
cultanes.inc <- read_dta("results/cultA_inc5_T_ALT.dta")
cultanes.edu <- read_dta("results/cultA_edu5_T_ALT.dta")


# Create DV indicator.
cultgss.egp$dv <- "Culture policy (GSS)"
cultanes.egp$dv <- "Culture policy (ANES)"

# Combine estimates for plotting.
# First need subset of EGP results.
egplst2 <- list(cultgss.egp, cultanes.egp)
egplst2 <- egplst2 %>%
  map(~ filter(., egp5==1 | egp5==2)) %>%
  map(~ rename(., b = b_egp)) %>%
  map(~ rename(., b_hi = b_egp_hi)) %>%
  map(~ rename(., b_lo = b_egp_lo)) %>%
  map(~ mutate(., var = ifelse(egp5 == 1, "Upper Class", 
                               "Upper Mid. Class")))

# Combine list of data frames.
time.cult <- bind_rows(egplst2)

# Create DV indicators and combine other data frames.
cultgsslst <- list(cultgss.subclass, cultgss.inc, cultgss.edu)
cultgsslst <- cultgsslst %>%
  map(~ mutate(., dv = "Culture policy (GSS)"))

cultaneslst <- list(cultanes.subclass, cultanes.inc, cultanes.edu)
cultaneslst <- cultaneslst %>%
  map(~ mutate(., dv = "Culture policy (ANES)"))

time.cult <- bind_rows(time.cult, cultgsslst, cultaneslst)

# Rename variables.
time.cult <- 
  time.cult %>% 
  mutate(var = recode_factor(var, 
                             "Upper Class" = "Upper Class",
                             "Upper Mid. Class" = "Upper Mid. Class",
                             "subclass" = "Subjective Class",
                             "inc5" = "Income",
                             "edu5" = "Education", .ordered = TRUE))

# Reorder DV.
dv.ord <- c("Culture policy (GSS)",
            "Culture policy (ANES)")
time.cult$dv <- factor(time.cult$dv, 
                       levels = dv.ord)


# Plots.

figTimeCult <- ggplot(time.cult, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "coral", size = 2) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF5.pdf", width=6, height=6)
figTimeCult
dev.off()



# Over time effects.
# Race.

# Load data.
# GSS.
racegss.egp <- read_dta("results/raceG_egp5_T_wo_ALT.dta")
racegss.subclass <- read_dta("results/raceG_subclass_T_wo_ALT.dta")
racegss.inc <- read_dta("results/raceG_inc5_T_wo_ALT.dta")
racegss.edu <- read_dta("results/raceG_edu5_T_wo_ALT.dta")

# ANES.
raceanes.egp <- read_dta("results/raceA_egp5_T_wo_ALT.dta")
raceanes.subclass <- read_dta("results/raceA_subclass_T_wo_ALT.dta")
raceanes.inc <- read_dta("results/raceA_inc5_T_wo_ALT.dta")
raceanes.edu <- read_dta("results/raceA_edu5_T_wo_ALT.dta")


# Create DV indicator.
racegss.egp$dv <- "Race policy (GSS)"
raceanes.egp$dv <- "Race policy (ANES)"

# Combine estimates for plotting.
# First need subset of EGP results.
egplst3 <- list(racegss.egp, raceanes.egp)
egplst3 <- egplst3 %>%
  map(~ filter(., egp5==1 | egp5==2)) %>%
  map(~ rename(., b = b_egp)) %>%
  map(~ rename(., b_hi = b_egp_hi)) %>%
  map(~ rename(., b_lo = b_egp_lo)) %>%
  map(~ mutate(., var = ifelse(egp5 == 1, "Upper Class", 
                               "Upper Mid. Class")))

# Combine list of data frames.
time.race <- bind_rows(egplst3)

# Create DV indicators and combine other data frames.
racegsslst <- list(racegss.subclass, racegss.inc, racegss.edu)
racegsslst <- racegsslst %>%
  map(~ mutate(., dv = "Race policy (GSS)"))

raceaneslst <- list(raceanes.subclass, raceanes.inc, raceanes.edu)
raceaneslst <- raceaneslst %>%
  map(~ mutate(., dv = "Race policy (ANES)"))

time.race <- bind_rows(time.race, racegsslst, raceaneslst)

# Rename variables.
time.race <- 
  time.race %>% 
  mutate(var = recode_factor(var, 
                             "Upper Class" = "Upper Class",
                             "Upper Mid. Class" = "Upper Mid. Class",
                             "subclass" = "Subjective Class",
                             "inc5" = "Income",
                             "edu5" = "Education", .ordered = TRUE))

# Reorder DV.
dv.ord <- c("Race policy (GSS)",
            "Race policy (ANES)")
time.race$dv <- factor(time.race$dv, 
                       levels = dv.ord)


# Plots.

figTimeRace <- ggplot(time.race, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "plum", size = 2) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF6.pdf", width=6, height=6)
figTimeRace
dev.off()




##################################################
# Estimates from separate regression models and
# random coefficient models plotted together.
##################################################

# Load MLM results.
# These results are created in file '20_figures_class.R'
vcls.mlm <- readRDS("results/time_vcls.rds")
redist.mlm <- readRDS("results/time_redist.rds")
cult.mlm <- readRDS("results/time_cult.rds")
race.mlm <- readRDS("results/time_race.rds")
vpol.mlm <- readRDS("results/time_vpol.rds")


# Separate regression model results.
# The 'rename_all' code removes underscores from variable names.
coeffs.gss <- read_dta("./results/coef_gss_sepmods.dta") %>% 
  rename_all(~ gsub("_", "", .))
coeffs.anes <- read_dta("./results/coef_anes_sepmods.dta") %>% 
  rename_all(~ gsub("_", "", .))

# Get subsets of data.
coeffs.gss <- filter(coeffs.gss, n >= 500)
coeffs.anes <- filter(coeffs.anes, n >= 500)

vcls.sep.gss <- filter(coeffs.gss, dv == "Vote Dem. (GSS)" 
                       & var != "Reduce Inequality"
                       & var != "Culture Policy"
                       & var != "Race Policy")
vcls.sep.anes <- filter(coeffs.anes, dv == "Vote Dem. (ANES)" 
                        & var != "Guarantee Jobs/Inc."
                        & var != "Culture Policy"
                        & var != "Race Policy")

vpol.sep.gss <- filter(coeffs.gss, var == "Reduce Inequality"
                       | var == "Culture Policy"
                       | var == "Race Policy")
vpol.sep.anes <- filter(coeffs.anes, var == "Guarantee Jobs/Inc."
                        | var == "Culture Policy"
                        | var == "Race Policy")

redist.sep.gss <- filter(coeffs.gss, dv == "Reduce inequality (GSS)")
redist.sep.anes <- filter(coeffs.anes, dv == "Guarantee jobs/inc. (ANES)")

cult.sep.gss <- filter(coeffs.gss, dv == "Culture policy (GSS)")
cult.sep.anes <- filter(coeffs.anes, dv == "Culture policy (ANES)")

race.sep.gss <- filter(coeffs.gss, dv == "Race policy (GSS)")
race.sep.anes <- filter(coeffs.anes, dv == "Race policy (ANES)")

# Combine separate model results.
vcls.sep <- bind_rows(vcls.sep.gss, vcls.sep.anes)
redist.sep <- bind_rows(redist.sep.gss, redist.sep.anes)
cult.sep <- bind_rows(cult.sep.gss, cult.sep.anes)
race.sep <- bind_rows(race.sep.gss, race.sep.anes)
vpol.sep <- bind_rows(vpol.sep.gss, vpol.sep.anes)


# Add new indicator.
vcls.mlm$model <- "Random coef. model"
redist.mlm$model <- "Random coef. model"
cult.mlm$model <- "Random coef. model"
race.mlm$model <- "Random coef. model"
vpol.mlm$model <- "Random coef. model"

# Rename and order factors.
vcls.sep <- vcls.sep %>%
  rename(b_lo = ll) %>%
  rename(b_hi = ul)

redist.sep <- redist.sep %>%
  rename(b_lo = ll) %>%
  rename(b_hi = ul)

cult.sep <- cult.sep %>%
  rename(b_lo = ll) %>%
  rename(b_hi = ul)

race.sep <- race.sep %>%
  rename(b_lo = ll) %>%
  rename(b_hi = ul)

vpol.sep <- vpol.sep %>%
  mutate(var2 = recode_factor(var,
                              "Reduce Inequality" = "Redist. pol.",
                              "Guarantee Jobs/Inc." = "Redist. pol.",
                              "Culture Policy" = "Culture pol.",
                              "Race Policy" = "Race pol.",
                              .ordered = TRUE)) %>%
  rename(b_lo = ll) %>%
  rename(b_hi = ul)


# Combine results.
time.vcls <- bind_rows(vcls.mlm, vcls.sep)
time.redist <- bind_rows(redist.mlm, redist.sep)
time.cult <- bind_rows(cult.mlm, cult.sep)
time.race <- bind_rows(race.mlm, race.sep)
time.vpol <- bind_rows(vpol.mlm, vpol.sep)

# Reorder factors.
time.vcls <- time.vcls %>%
  mutate(dv = fct_relevel(dv,
                          "Vote Dem. (GSS)",
                          "Vote Dem. (ANES)")) %>%
  mutate(var = fct_relevel(var,
                           "Upper Class",
                           "Upper Mid. Class",
                           "Subjective Class",
                           "Income",
                           "Education"))

time.redist <- time.redist %>%
  mutate(dv = fct_relevel(dv,
                          "Reduce inequality (GSS)",
                          "Guarantee jobs/inc. (ANES)")) %>%
  mutate(var = fct_relevel(var,
                           "Upper Class",
                           "Upper Mid. Class",
                           "Subjective Class",
                           "Income",
                           "Education"))

time.cult <- time.cult %>%
  mutate(dv = fct_relevel(dv,
                          "Culture policy (GSS)",
                          "Culture policy (ANES)")) %>%
  mutate(var = fct_relevel(var,
                           "Upper Class",
                           "Upper Mid. Class",
                           "Subjective Class",
                           "Income",
                           "Education"))

time.race <- time.race %>%
  mutate(dv = fct_relevel(dv,
                          "Race policy (GSS)",
                          "Race policy (ANES)")) %>%
  mutate(var = fct_relevel(var,
                           "Upper Class",
                           "Upper Mid. Class",
                           "Subjective Class",
                           "Income",
                           "Education"))


# Plots.

# Over time effects.
# Voting and class.

figTimeVcls <- ggplot(time.vcls, aes(y = b, x = year, color = model, shape = model)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 0.75, color = "gray", width=0,
                position = position_dodge(width = 1.5)) +
  geom_point(size = 1.5,
             position = position_dodge(width = 1.5)) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients",
       color = "", shape = "") +
  fte_theme() +
  theme(legend.position = "bottom") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF7.pdf", width=6, height=6)
figTimeVcls
dev.off()


# Over time effects.
# Redistribution.

figTimeRedist <- ggplot(time.redist, aes(y = b, x = year, color = model, shape = model)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 0.75, color = "gray", width=0,
                position = position_dodge(width = 1.5)) +
  geom_point(size = 1.5,
             position = position_dodge(width = 1.5)) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients",
       color = "", shape = "") +
  fte_theme() +
  theme(legend.position = "bottom") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white")) 


pdf(file = "figF9.pdf", width=6, height=6)
figTimeRedist
dev.off()


# Over time effects.
# Culture policy.

figTimeCulture <- ggplot(time.cult, aes(y = b, x = year, color = model, shape = model)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 0.75, color = "gray", width=0,
                position = position_dodge(width = 1.5)) +
  geom_point(size = 1.5,
             position = position_dodge(width = 1.5)) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients",
       color = "", shape = "") +
  fte_theme() +
  theme(legend.position = "bottom") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white")) 


pdf(file = "figF10.pdf", width=6, height=6)
figTimeCulture
dev.off()


# Over time effects.
# Race policy.

figTimeRace <- ggplot(time.race, aes(y = b, x = year, color = model, shape = model)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 0.75, color = "gray", width=0,
                position = position_dodge(width = 1.5)) +
  geom_point(size = 1.5,
             position = position_dodge(width = 1.5)) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients",
       color = "", shape = "") +
  fte_theme() +
  theme(legend.position = "bottom") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white")) 


pdf(file = "figF11.pdf", width=6, height=6)
figTimeRace
dev.off()


# Over time effects.
# Voting and policy.

figTimeVpol <- ggplot(time.vpol, aes(y = b, x = year, color = model, shape = model)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 0.75, color = "gray", width=0,
                position = position_dodge(width = 1.5)) +
  geom_point(size = 1.5,
             position = position_dodge(width = 1.5)) +
  facet_grid(var2 ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients",
       color = "", shape = "") +
  fte_theme() +
  theme(legend.position = "bottom") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white")) 


pdf(file = "figF14.pdf", width=6, height=6)
figTimeVpol
dev.off()



##################################################
# White only subsample models.
##################################################

# Voting over time.

# Load data.
# GSS.
vgss.egp.wo <- read_dta("./results/dvote_egp1_gss_time_wo.dta")
vgss.subclass.wo <- read_dta("./results/dvote_subclass_gss_time_wo.dta")
vgss.inc.wo <- read_dta("./results/dvote_inc5_gss_time_wo.dta")
vgss.edu.wo <- read_dta("./results/dvote_edu5_gss_time_wo.dta")

# ANES
vanes.egp.wo <- read_dta("./results/dvote_egp1_anes_time_wo.dta")
vanes.subclass.wo <- read_dta("./results/dvote_subclass_anes_time_wo.dta")
vanes.inc.wo <- read_dta("./results/dvote_inc5_anes_time_wo.dta")
vanes.edu.wo <- read_dta("./results/dvote_edu5_anes_time_wo.dta")

# Create DV indicator.
vgss.egp.wo$dv <- "Vote Dem. (GSS)"
vanes.egp.wo$dv <- "Vote Dem. (ANES)"

# Combine estimates for plotting.
# First need to rename and create variables.
vegplst.wo <- list(vgss.egp.wo, vanes.egp.wo)
vegplst.wo <- vegplst.wo %>%
  map(~ rename(., b = b_egp)) %>%
  map(~ rename(., b_hi = b_egp_hi)) %>%
  map(~ rename(., b_lo = b_egp_lo)) %>%
  map(~ mutate(., var = ifelse(egp5 == 1, "Upper Class", 
                               "Upper Mid. Class")))

# Combine list of data frames.
time.vote.wo <- bind_rows(vegplst.wo)

# Create DV indicators and combine other data frames.
vglst.wo <- list(vgss.subclass.wo, vgss.inc.wo, vgss.edu.wo)
vglst.wo <- vglst.wo %>%
  map(~ mutate(., dv = "Vote Dem. (GSS)"))

valst.wo <- list(vanes.subclass.wo, vanes.inc.wo, vanes.edu.wo)
valst.wo <- valst.wo %>%
  map(~ mutate(., dv = "Vote Dem. (ANES)"))

time.vote.wo <- bind_rows(time.vote.wo, vglst.wo, valst.wo)

# Rename variables.
time.vote.wo <- 
  time.vote.wo %>% 
  mutate(var = recode_factor(var, 
                             "Upper Class" = "Upper Class",
                             "Upper Mid. Class" = "Upper Mid. Class",
                             "subclass" = "Subjective Class",
                             "inc5" = "Income",
                             "edu5" = "Education", .ordered = TRUE))

# Reorder DV.
vdv.ord <- c("Vote Dem. (GSS)",
             "Vote Dem. (ANES)")
time.vote.wo$dv <- factor(time.vote.wo$dv, 
                          levels = vdv.ord)


# Plots.

figTimeVoteWO <- ggplot(time.vote.wo, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "#00ba38", size = 2) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF8.pdf", width=6, height=6)
figTimeVoteWO
dev.off()


# Redistribution over time.

# Load data.
# Inequality.
ineq.egp.wo <- read_dta("./results/ineq_egp5_time_wo.dta")
ineq.subclass.wo <- read_dta("./results/ineq_subclass_time_wo.dta")
ineq.inc.wo <- read_dta("./results/ineq_inc5_time_wo.dta")
ineq.edu.wo <- read_dta("./results/ineq_edu5_time_wo.dta")

# Jobs and income.
jobs.egp.wo <- read_dta("./results/jobs_egp5_time_wo.dta")
jobs.subclass.wo <- read_dta("./results/jobs_subclass_time_wo.dta")
jobs.inc.wo <- read_dta("./results/jobs_inc5_time_wo.dta")
jobs.edu.wo <- read_dta("./results/jobs_edu5_time_wo.dta")


# Create DV indicator.
ineq.egp.wo$dv <- "Reduce inequality (GSS)"
jobs.egp.wo$dv <- "Guarantee jobs/inc. (ANES)"

# Combine estimates for plotting.
# First need subset of EGP results.
egplst.wo <- list(ineq.egp.wo, jobs.egp.wo)
egplst.wo <- egplst.wo %>%
  map(~ filter(., egp5==1 | egp5==2)) %>%
  map(~ rename(., b = b_egp)) %>%
  map(~ rename(., b_hi = b_egp_hi)) %>%
  map(~ rename(., b_lo = b_egp_lo)) %>%
  map(~ mutate(., var = ifelse(egp5 == 1, "Upper Class", 
                               "Upper Mid. Class")))

# Combine list of data frames.
time.redist.wo <- bind_rows(egplst.wo)

# Create DV indicators and combine other data frames.
ineqlst.wo <- list(ineq.subclass.wo, ineq.inc.wo, ineq.edu.wo)
ineqlst.wo <- ineqlst.wo %>%
  map(~ mutate(., dv = "Reduce inequality (GSS)"))

jobslst.wo <- list(jobs.subclass.wo, jobs.inc.wo, jobs.edu.wo)
jobslst.wo <- jobslst.wo %>%
  map(~ mutate(., dv = "Guarantee jobs/inc. (ANES)"))

time.redist.wo <- bind_rows(time.redist.wo, ineqlst.wo, jobslst.wo)

# Rename variables.
time.redist.wo <- 
  time.redist.wo %>% 
  mutate(var = recode_factor(var, 
                             "Upper Class" = "Upper Class",
                             "Upper Mid. Class" = "Upper Mid. Class",
                             "subclass" = "Subjective Class",
                             "inc5" = "Income",
                             "edu5" = "Education", .ordered = TRUE))

# Reorder DV.
dv.ord <- c("Reduce inequality (GSS)",
            "Guarantee jobs/inc. (ANES)")
time.redist.wo$dv <- factor(time.redist.wo$dv, 
                            levels = dv.ord)


# Plots.

figTimeRedistWO <- ggplot(time.redist.wo, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "#619cff", size = 2) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF12.pdf", width=6, height=6)
figTimeRedistWO
dev.off()


# Culture policy over time.

# Load data.
# GSS culture.
cultgss.egp.wo <- read_dta("./results/cultgss_egp5_time_wo.dta")
cultgss.subclass.wo <- read_dta("./results/cultgss_subclass_time_wo.dta")
cultgss.inc.wo <- read_dta("./results/cultgss_inc5_time_wo.dta")
cultgss.edu.wo <- read_dta("./results/cultgss_edu5_time_wo.dta")

# ANES culture.
cultanes.egp.wo <- read_dta("./results/cultanes_egp5_time_wo.dta")
cultanes.subclass.wo <- read_dta("./results/cultanes_subclass_time_wo.dta")
cultanes.inc.wo <- read_dta("./results/cultanes_inc5_time_wo.dta")
cultanes.edu.wo <- read_dta("./results/cultanes_edu5_time_wo.dta")


# Create DV indicator.
cultgss.egp.wo$dv <- "Culture policy (GSS)"
cultanes.egp.wo$dv <- "Culture policy (ANES)"

# Combine estimates for plotting.
# First need subset of EGP results.
egplst2.wo <- list(cultgss.egp.wo, cultanes.egp.wo)
egplst2.wo <- egplst2.wo %>%
  map(~ filter(., egp5==1 | egp5==2)) %>%
  map(~ rename(., b = b_egp)) %>%
  map(~ rename(., b_hi = b_egp_hi)) %>%
  map(~ rename(., b_lo = b_egp_lo)) %>%
  map(~ mutate(., var = ifelse(egp5 == 1, "Upper Class", 
                               "Upper Mid. Class")))

# Combine list of data frames.
time.cult.wo <- bind_rows(egplst2.wo)

# Create DV indicators and combine other data frames.
cultgsslst.wo <- list(cultgss.subclass.wo, cultgss.inc.wo, cultgss.edu.wo)
cultgsslst.wo <- cultgsslst.wo %>%
  map(~ mutate(., dv = "Culture policy (GSS)"))

cultaneslst.wo <- list(cultanes.subclass.wo, cultanes.inc.wo, cultanes.edu.wo)
cultaneslst.wo <- cultaneslst.wo %>%
  map(~ mutate(., dv = "Culture policy (ANES)"))

time.cult.wo <- bind_rows(time.cult.wo, cultgsslst.wo, cultaneslst.wo)

# Rename variables.
time.cult.wo <- time.cult.wo %>% 
  mutate(var = recode_factor(var, 
                             "Upper Class" = "Upper Class",
                             "Upper Mid. Class" = "Upper Mid. Class",
                             "subclass" = "Subjective Class",
                             "inc5" = "Income",
                             "edu5" = "Education", .ordered = TRUE))

# Reorder DV.
dv.ord <- c("Culture policy (GSS)",
            "Culture policy (ANES)")
time.cult.wo$dv <- factor(time.cult.wo$dv, 
                          levels = dv.ord)


# Plots.

figTimeCultWO <- ggplot(time.cult.wo, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "coral", size = 2) +
  facet_grid(var ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF13.pdf", width=6, height=6)
figTimeCultWO
dev.off()



# Voting and policy over time.

# Load data.

dvptime.gss.wo <- read_dta("./results/dvote_policy_gss_time_wo.dta")
dvptime.anes.wo <- read_dta("./results/dvote_policy_anes_time_wo.dta")

# Create DV indicator.
dvptime.gss.wo$dv <- "Vote Dem. (GSS)"
dvptime.anes.wo$dv <- "Vote Dem. (ANES)"

# Combine results.
dvptime.all.wo <- bind_rows(dvptime.gss.wo, dvptime.anes.wo)

# New label var.
dvptime.all.wo <- dvptime.all.wo %>% 
  mutate(var2 = recode_factor(var, 
                              "govineq01"="Redist. pol.",
                              "jobs01"="Redist. pol.",
                              "culpol01"="Culture pol.",
                              "racepol01"="Race pol.", 
                              .ordered = TRUE))

# Reorder factors.
dvptime.all.wo$var2 <- factor(dvptime.all.wo$var2, 
                              levels = c("Redist. pol.", 
                                         "Culture pol.",
                                         "Race pol."))
dvptime.all.wo$dv <- factor(dvptime.all.wo$dv, 
                            levels = c("Vote Dem. (GSS)",
                                       "Vote Dem. (ANES)"))


# Plots.

figTimeDVPwo <- ggplot(dvptime.all.wo, aes(y = b, x = year)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_errorbar(aes(ymin = b_lo, ymax = b_hi), 
                size = 1, color = "gray", width=0) +
  geom_point(color = "palegreen3", size = 2) +
  facet_grid(var2 ~ dv, scales = "free_y") +
  labs(x = "Year", 
       y = "Estimated Coefficients") +
  fte_theme() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1972, 2018, 8)) +
  theme(plot.background = element_rect(fill = "white"))

pdf(file = "figF15.pdf", width=6, height=6)
figTimeDVPwo
dev.off()





